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Array deconvolution is commonly used in aeroacoustic analysis to remove the influ- 
ence of a microphone array’s point spread function from a conventional beamforming map. 
Unfortunately, the majority of deconvolution algorithms assume that the acoustic sources 
in a measurement are incoherent, which can be problematic for some aeroacoustic phe- 
nomena with coherent, spatially-distributed characteristics. While several algorithms have 
been proposed to handle coherent sources, some are computationally intractable for many 
problems while others require restrictive assumptions about the source field. Newer gen- 
eralized inverse techniques hold promise, but are still under investigation for general use. 
An alternate coherent deconvolution method is proposed based on a wavespace transfor- 
mation of the array data. Wavespace analysis offers advantages over curved-wave array 
processing, such as providing an explicit shift-invariance in the convolution of the array 
sampling function with the acoustic wave field. However, usage of the wavespace transfor- 
mation assumes the acoustic wave field is accurately approximated as a superposition of 
plane wave fields, regardless of true wavefront curvature. The wavespace technique lever- 
ages Fourier transforms to quickly evaluate a shift- invariant convolution. The method is 
derived for and applied to ideal incoherent and coherent plane wave fields to demonstrate 
its ability to determine magnitude and relative phase of multiple coherent sources. Multi- 
scale processing is explored as a means of accelerating solution convergence. A case with a 
spherical wave front is evaluated. Finally, a trailing edge noise experiment case is consid- 
ered. Results show the method successfully deconvolves incoherent, partially-coherent, and 
coherent plane wave fields to a degree necessary for quantitative evaluation. Curved wave 
front cases warrant further investigation. A potential extension to nearfield beamforming 
is proposed. 


Nomenclature 


Dimensions 


kx i foy 
, k^ 


%,y 


Wavespace dimensions 
Conjugate-power wavespace dimensions 
Spatial dimensions 
Conjugate-power spatial dimensions 


Symbols 

Relaxation parameter 
Example system matrix 
Matrix used in iteration stability analysis 
Matrix used in iteration stability analysis 
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Acronym for Computational Fluid Dynamics 

Acronym for Cross-Spectral Matrix 

Isentropic speed of sound 

Matrix used in iteration stability analysis 

Vector of wavespace transform terms 

Expected value operation 

Normalized solution error 

Temporal frequency 

Acronym for Fast Fourier Transform 

Example cross-spectral matrix 

Acronym for Graphics Processor Unit 

Identity matrix 

imaginary unit 

Maximum extent of wavespace grid 

Total number of microphones in array 

Number of grid points in a beam map 

Number of grid points in a given dimension in wavespace 

Number of grid points in a given dimension required for linear convolution 

Sampled pressure field on array face 

Pressure sampled at microphone m for a given narrowband frequency 

Wavespace transform of sampled pressure field 

Spatial cross-spectrum 

Measured wavespace cross-spectrum 

Array point spread function 

Continuous pressure field on array face 

Wavespace transform of true pressure field 

Wavespace cross-spectrum of true pressure field 

Forward transform of Q 

Spectral radius of a matrix 

Solution residual 

Convolution of source field with array response 

Forward transform of R 

Array sampling function in space 

Wavespace transform of array sampling function 

Wavespace cross-power sampling function 

Forward transform of S 

Vector of example sources 

Vector of example observations 

Grid point spacing in wavespace 

Dirac delta function 

Wavespace averaging function 


Subscripts and Superscripts 

() u Index of average bin in wavespace 

() H Hermitian transpose 

()« Current iteration number 

()* Complex conjugate 

()' Shifted coordinate 

() Wavespace average 
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I. Introduction 


T he limitations of conventional beamforming impede the ability to extract quantitative information from 
traditional beam maps. Simple integration techniques can be used to provide approximate field values, 
but in general the inverse array problem of deconvolution must be solved to extract quantitative information 
from array data. Many frequency- domain deconvolution algorithms exist for incoherent source fields, for ex- 
ample DAMAS,® DAMAS2plSC-DAMASP and CLEAN-SC. 3 4 ^ However, when a source field contains regions 
of coherence, there are fewer algorithm selections. Methods such as MACS are computationally reasonable, 
but assume a sparse source field. 5 Similarly, LORE is a quick algorithm but only solves for a subset of the 
scan grid. 0 ] DAMAS-C has the potential to evaluate coherent sources over the entire scan grid of interest, but 
to do so is computationally challenging, even for small source fields. 7 Generalized inverse methods involving 
either eigenvalue subsets of the CS A/® or the entire CSM of an array data setP may resolve distributed, 
partially-coherent source fields, but the range of applicability and limitations of these methods have yet to 
be completely addressed. 

Deconvolution is fundamentally the inversion of a linear system, expressed in Eq. 0 , with the application 
of constraints. For a vector of observations y and a matrix A modeling acoustic propagation, assuming the 
observations are uncontaminated by measurement noise, deconvolution attempts to solve 

y = Ax (1) 

for cc, representing the sources which pass through the system matrix to generate the observations. For 
many problems of interest, A is nearly singular, so direct inversion may not be possible, or yield a false 
solution. Algorithms such as DAMAS and DAMAS2 handle this inversion process iteratively, where an 
initial guess of x is multiplied by A and then subtracted from y to compute a residual. This residual is used 
to update the guess of x. DAMAS2 makes a shift-invariant assumption, allowing the use of a Fourier-based 
technique to perform a fast convolution of x with A? Individual iterations of DAMAS2 are significantly 
faster than DAMAS for a given problem, but depending on the validity of the shift- invariant assumption and 
desired convergence criteria, more iterations of DAMAS2 may be required and accuracy may be limited. A 
smoothing filter may be used with DAMAS 2 to favor certain solution sets and accelerate convergence. 

Coherent deconvolution, where elements of x are statistically related to each other, is also modeled by 
Eq. 0 . In this case, cross-spectral terms between grid points are included in x and i/, and A updated 
accordingly. Constraints on these additional elements are defined by a cross-spectral inequalityp 0 rather 
than a positivity constraint, and in general the terms will be complex. Unfortunately, scaling of general 
coherent deconvolution becomes problematic. If a beam map of interest has n g grid points, x and y each 
have n 2 g elements, and A has n 4 elements. For example, a 20 x 20 beam map has n g = 400 grid points. This 
means x and y each have 20 4 elements to account for all possible auto- and cross-spectral terms, meaning 
A has 20 8 , or 25.6 billion, elements. Conjugate- symmetry reduces the number of independent cross-spectral 
terms by a factor of two, but the scaling remains. If cross-spectral relationships between all of the grid points 
are considered, A is generally not sparse. While such a matrix of this size is possible to generate and store, 
it does not lend itself to efficient computation nor to more realistic grids. 

To make the evaluation of Ax more tractable for this form of deconvolution, it would be helpful to 
define the problem such that every element of A does not need to be computed and stored. Additionally, 
a fast convolution technique similar to that used by DAMAS2 is desired. A formulation of the coherent 
deconvolution process which does this is presented. A shift- invariance assumption, where the element of A 
relating a given element of y to an element of x is dependent only on the spatial separation of the two points, 
is applied to the coherence problem, just as with DAMAS 2. This is done by analyzing the array data in 
wavenumber space, equivalent to plane wave beamforming. The analysis sacrifices the depth perception of 
nearfield beamforming due to the plane wave assumption. 

The general problem of a transformation of a coherent pressure field to wavespace is first discussed in 
Section II. Basic implementation details are then addressed in Section III. Potential methods for accelerating 
convergence through faster computation or iteration count reduction are considered in Section IV. Section 
V presents results for different combinations of incoherent and coherent plane wave fields. Section VI 
evaluates a simulated curved wave field generated by a point source, as well as data from a trailing edge noise 
measurement. Finally, a method for nearfield beamforming involving the transformation of the deconvolved 
wavespace spectrum is proposed. 
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II. Formulation 


A formulation of the wavespace deconvolution problem which follows the form of Eq. 0 is necessary. 
Here, y is the set of array observations transformed to wavespace, and A is the model of the array sam- 
pling function. This array sampling function in wavespace is first constructed in regular space in terms of 
sampling theory, such that the true pressure field in space is multiplied by a distribution of delta functions 
corresponding to microphone locations in a planar array? 11 This multiplication in the spatial domain be- 
comes an equivalent convolution in the Fourier transform of the spatial domain, or the wavespace domain. 
The sampling function operating on the pressure field on a planar array face is given by 

M 

s(x,y) = E S(x-x m ,y-y m ). (2) 

m= 1 

The sampled, temporally-Fourier transformed narrowband pressure field (frequency-dependence suppressed) 
is then given by 


p(x,y) = s (x, y) q (x, y) 

M 

^ ^ Pmfi %rn i V Vm) • 

m= 1 

The 2-D spatial Fourier transforms of these quantities are thus 

r r i ^ 

s(k x ,k y ) = s (x,y) e~ j ^ kxX+kyV ^dxdy = — ^ e -o(k x xm+ky yrn ) 

' ' m= 1 


(3) 


(4) 


and 


P (k x , ky ) = JJ S (k x -k’ x ,ky- ky) q (k' x , k' y ) dk' x dk' y 

1 M 

— _ V r> p-j( k xXm+k y ym ) 

- M 2^ Pme 


(5) 


m= 1 


where the normalization by M is included such that a unit-magnitude plane wave in physical space has unit 
magnitude in wave space. Note that no k z transform is performed as all microphones are assumed co-planar. 

Eq. 0 is a formulation of the basic wavespace sampling problem for a single block of data. To handle 
the stochastic nature of the pressure field, auto- and cross-powers are used to represent the ensemble-average 
characteristics of the field. A CSM contains the auto- and cross-powers of the array microphones in the 
frequency domain ensemble- averaged across many blocks of data, and is represented by 


G = 


E \p* (x 1 ,y l )p{x 1 ,y 1 )} 
E\p* (x2,y 2 )p(xi,yi)] 


E \p* (x 1 ,y 1 )p(x 2 ,y 2 )] 
E \p* (x 2 ,y 2 )p(x 2 ,y 2 )} 


E\p* (xi ,yi)p(x M ,y M )} 
E\p* ( x 2 ,y 2 )p(x M ,y M )] 


(6) 


_E \p* (x M ,yM)p(xi,yi)} E \p* (xm ,Pm) P {x 2 ,y 2 )\ ... E\p* (x M ,yM)p(xM,yM)\_ 

Note that the 2-D CSM contains 4-D spatial covariance information, 

P (£, V, x,y)=E \p* (^, rf) p (x, y)\ . 

The CSM could be directly transformed to wavespace through a 4-D Fourier transform, 

$ { P (£, r/,x,y)j = JJJJp (£, V, x , y) e ~^ k ^ +k ^ +k - x+k « y '>d^dr ] dxdy. 

However, for this discussion the wavespace cross-spectrum is the desired quantity, which is defined by 

P (k^, ky) k x , ky) — E \p (&£, ky) P {k x , ky )] . 


(7) 


(8) 


(9) 
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This can be determined by substituting the first line of Eq. (J5| into Eq. ([9|, 

P (&£, ky, k x , ky) = E [p (k^, ky) p ( k x , /c?/)] 


= E 


JJ s* (fc c - k v - kQ q* (*£, kQ dkpk'yX 

JJ s (k x — k' x , k y — ky) q (k' x , k' y ) dk' x dk' y (10) 

= JJJJ s * - k v - K) s ( kx - k '*’ k v - K ) x 

E [r (fy, k'y) q (k' x , k'y)] 

= JUJ S (fc € - ky ~ k'y , k X ~ K, ky ~ k'y) Q (%, k'y , ^ . 

This shows that the measured wavespace cross-spectrum is the convolution of the wavespace cross-power 
sampling function S and the true wavespace cross-spectrum, Q. 

The measured wavespace cross-spectrum P can be computed by substituting the second line of Eq. © 
into Eq. ([ 9 ]), giving 


P (&£ , ky, k x , ky) — E [p (/c^ , ky) p ( k x , ky)] 


= E 


' M 1 M 

_jL n* JfaZn+knrin) 

M M ^ 


Pm? 


- j{k x Xm+k y ym ) 


n=l 
M M 




m= 1 


(fc££n + ^^n) p -i(fc®®m+fc y 2/m) 


( 11 ) 


n=l m=l 


This is similar to the derivation given by CaporlEI which relates the cross-spectrum between two points 
to their separation in physical space and the wavespace spectrum, but in a form of a forward transform 
instead of an inverse. As a 4D representation of 2D wavespace, this representation captures the coherence 
relationship between wavespace locations. The term E [p^Pm] in Eq. is an entry in the CSM. Thus, 
the double-summation can be re-expressed as 


P(k i ,ky,k x ,ky) = T e H Ge, 


with the Hermitian transpose of e given by 

e j(k^2-\-k v rj 2 ) 


e H = 


pj ~\~krjT] M ) 


( 12 ) 


(13) 


Similarly, the wavespace cross-power sampling function can be constructed as 


S(k^ky,k X ,ky) ^ 


1 1 
1 1 

1 1 


e. 


(14) 


Note that when diagonal removal is applied, the normalization should change from M 2 to M 2 — M to 
maintain unit response in the transformation as m = n terms are neglected. 

As previously stated, formulation of the problem is desired which follows Eq. ©• In this section, a 
model is constructed in which P contains the information in the observation vector y, and S contains the 
information in the system matrix A. Taking the source field Q as the solution vector x , the problem of 
interest is to find the source field which best fits the data in the observation vector, obeying the given 
constraints and avoiding explicit computation or storage of the entirety of A. 

III. Implementation 

To solve this problem, P and S are first constructed for a discrete set of coordinates in wavespace. 
For simplicity, a hypercube grid is considered, with each coordinate set spanning k^ = ky = k x = k y = 
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[—k max : Afc : k max \ for a total grid size of rik x rik x x n^. The selection of the parameters k max , Afc 
and rik are discussed subsequently, but are dependent on the array, the problem of interest and available 
computational resources. 

An iterative solution of Eq. 0 can take the form of 

X (i+1) = x (i) + I (y - AbW) . (15) 

Constraints are enforced after every iteration. The general constraints considered for DAMAS-C are adapted. 6 7 
A real positivity constraint is applied to autospectral terms, 


5R {<2 (k x , k y , k x , k y )| >0 (16) 

^ (J^ X , ky^ k X: ky^) J* — 0. 

Cross- spectral terms must obey the cross-spectral inequality and follow appropriate conjugate behavior, 


Q (&£, ky, k x , ky ) ^ Q ( k x 5 ky , A^, A)y) Q (Aj^, ky, k £, A)^) 

Q ( k x i ky , A)£ , A^) — (2 (&£ , ky , k x , A)^) . 


(17) 


After each iteration, autospectral constraints are applied first, followed by cross-spectral constraints. 

The algorithm takes the same structure as that used with DAM AS 2, where the equivalent evaluation 
of Ax^ from Eq. (15) is computed as a function evaluation through the use of Fourier transforms, rather 
than an explicit linear algebra operation.^ Unlike DAMAS2, no regularization filter is applied at this stage 
of research, as the effect this may have on phase relationships is unclear. 


1. Forward-transform S to S via a 4-D transform like Eq. 0- 

2. Initialize Q to zeroes. 


3. Iterate (to convergence or maximum desired function evaluations) 

• Forward transform QW to Q , again using Eq. 

• Compute the product R = SQ. 

• Inverse transform R to R to obtain the effective convolution of S with Q^\ 

• Update = QW + ± (p - 

• Enforce constraints. 


Note that for linear convolution, Q must be zero-padded, and S must be constructed on the larger, padded 
grid scale. Here, the grid coordinates all span [— 2k max : Ak : 2 k max \, for a total grid size of n s xn s xn s xn s . 
If rik is selected to be odd such that the plane wave bin Q (0, 0, 0, 0) is in the solution set, then n s = 2 rik — 1 
and the total grid size is (2 rik — 1) 4 ? or approximately 16 times the size of the baseline grid. Q must be 
padded to this size prior to convolution, and the appropriate subset of grid points selected in the P — R step 
of the algorithm. 

The relaxation parameter, a, is a critical component of the algorithm, as the iterative update is completely 
unstable without it. With DAMAS2, this parameter is specified as the sum of the absolute value of the 2D 
array response within the baseline grid domain. Extended to this 4D problem, it is computed by 


— ^A ^A S (A^J ky, k x , ky) 


/c^ — k max k r j — kmax kx — ^max ky — k^ 


(18) 


This value is effectively the maximum row sum of absolute values of the A representation of S in Eq. (15). 
This corresponds to the matrix norm HA^. 14 As the spectral radius of a square matrix must be less than or 
equal to its operator norms, normalizing the iterative process by this norm is a conservative way to stabilize 
the solution procedure. 
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While this algorithm bypasses explicit storage of A, which would have a size of n|, problem scaling can 
still be an issue. Specifically, n 4 problem scaling is still non-trivial. For example, if = 25, each n 4 -sized 
array, where n s =49, contains 5.8 million elements. While this has moderate storage requirements, 46 MB 
for single-precision complex, it is still a significant amount of data to handle with an algorithm which may 
require a significant number of iterations. Hardware-tuned FFT algorithms® can help mitigate this expense, 
but the cost of a given problem size is still a driving concern, since based on Eq. (18), larger problems will 
in general require greater relaxation parameters for a given array resolution. 

Previous research 1 1 indicates that in deconvolution the grid spacing should be between approximately 5% 
and 20% of the sampling pattern’s 3-dB main lobe. In wavespace, this sets an effective maximum value for 
A k for a given array. The maximum wavenumber k max of interest is highly problem-dependent. In a limiting 
case of plane waves arriving within a limited cone of directions, k max can simply be selected to properly 
encompass the appropriate region within wavespace. However, for a general problem where no information 
is known about the acoustic field k max will need to meet, if not exceed, the acoustic wavenumber, 2irf /cq 
(rad/m), or //co (m -1 ). These constraints, combined with the aforementioned scaling issue, currently limit 
the method to either problems with limited direction of arrival, or low frequencies. For example, the 0.74 m 
aperture outer array used at the University of Florida Aeroacoustic Flow Facility (UFAFFjl 6 can adequately 
capture the entire acoustic radiation circle with a A k of 17% of the sampling pattern’s 3-dB main lobe, 2 
m -1 , at / = 1 kHz with = 17. An individual transform of a grid required for this scale is reasonably 
quick, and thus can be considered for analysis. Higher frequencies become an issue, as acoustic wavenumber 
scales linearly with frequency. To maintain a given Afc, this means that must also scale linearly with 
frequency, so the overall problem scales as the fourth power of frequency. Similarly, extending A: maa , beyond 
the acoustic wavenumber to evaluate subsonic hydrodynamic contamination can only be considered for very 
low frequencies. 

As formulated above and tested, the algorithm worked for simple simulated problems. However, for a grid 
with rik = 9, approximately 500,000 iterations were required to converge for an ideal point source, defined 
by 

err (*+l) _ ||p- . (19) 

The residual change in solution is also followed, defined by, 

res (i +!) = || Qb +1 ) - QW||, (20) 

Numeric convergence is achieved when the error reaches the minimum unit-round based on floating-point sin- 
gle precision. 14 i For larger grids, numeric convergence never occurred, although reasonable- looking solutions 
were achieved. For these, the residual was tracked until it became small relative to the solution magnitude. 
Changing to double-precision did not appear to improve convergence behavior for the = 9 grids, and 
came with a significant run time and memory cost. At this stage of research, the effect of floating-point 
precision on the convergence behavior of larger grids was not considered. 

IV. Solution Acceleration 

To improve the utility of this deconvolution technique, acceleration methods must be explored. Three 
possibilities are considered: reducing the operating time of an individual iteration, using a more aggressive 
relaxation parameter a, and improving the initial solution guess Q(°). 


Computation Time 

A direct reduction in iteration time is easiest to consider, as 500,000 iterations may appear reasonable if the 
total code run time is 10 seconds. Hardware-optimized FFT algorithms are continuously improving, and 
some companies are now offering pre-made libraries for GPU computing^ 17 potentially providing significant 
performance improvements. Additionally, pruned FFT methods can be considered. The majority of the 
elements in Q are zeros to meet padding requirements. The majority of the elements in R are thrown away 
prior to subtraction from P. Hardware implementations of pruned FFTs existf® and may provide significant 
benefits in cases of 4-dimensional zero-padding. However, neither these nor GPU methods are explored at 
this stage of the research. 

One potential algorithm was considered. The 4-dimensional FFT is, as previously mentioned, a fast 
mechanism for evaluating the Ax^ term from Eq. (15) without storing the entirety of A. Alternative fast 
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convolution techniques exist. One such technique takes advantage of the fact that, with appropriate element 
ordering, A has a multilevel block- Toeplitz structure. It is possible to store only the unique elements of this 
matrix in a vector and perform a fast convolution with an appropriately-padded x , providing an alternative 
storage- free fast convolution. 19 The algorithm as cited was implemented, with some modifications to the 
zero-padding code and element access for performance. It was found to perform slightly better for some grid 
sizes and slightly worse for others. While it may hold some promise for improvement, the method is not 
further explored here. 


Relaxation 


The relaxation parameter computed in Eq. (18) is a conservative one. With some effort, a more aggressive 
parameter which maintains stability may be derived, following the method presented by Atkinson. 14 Eq. 
0 can be re-structured by splitting A, 

A = B-C, (21) 


and re-writing as 

The iterative method is thus given as 
Matrix D can be defined as 


Bx = y + Cx. 

Bx {i+1) = y + Cx «. 


D = B'C. 

The stability criterion for the (unconstrained) iteration method is then given as 

rv (D) < 1. 


( 22 ) 

(23) 

(24) 

(25) 


With some rearrangement, Eq. (15) yields B = al and C = al — A. This gives D = I — -A. Since 


r a (D) is the magnitude of the largest eigenvalue of D , finding this eigenvalue and iterating a until the 
spectral radius drops below unity should yield a more aggressive value for a. This may appear problematic 
at first as the goal of this method is to avoid computing and storing matrices of size A, which would be 
necessary to determine all of the eigenvalues of D. However, there are methods for computing the largest 
eigenvalue of D as a function call evaluating the matrix vector product Dx for a given input vector cc, rather 
than using D itself. Specifically, the MATLAB function eigs can output the largest eigenvalue of a large, 
dense, non-symmetric complex matrix modeled by a function call. As used with these function calls, eigs 
fails to converge when the eigenvalue of D drops below unity, so the smallest value of a for which eigs fails 
to give a converged output is used. For this study, the value varies between one-half and one-twentieth of 
the value given by Eq. (18) for the modeled arrays and grids, and appears stable for the cases considered. 


The reduced relaxation parameter significantly improves the convergence behavior of the algorithm. Cases 
with rik = 9 reach acceptable convergence in 50,000 iterations rather than 500,000. 

Initial Guess 

Additional improvement is desired for larger grids. An improved initial guess can have a significant impact 
on the direction and rate of convergence of the method, and is worth investigating. A variety of initialization 
methods were considered (e.g. Wiener filtering, initializing with DAMAS3P starting from P). However, 
these techniques yielded behavior generally worse than starting with Q — 0. However, one method which 
did show improvement for some situations is multiscale analysis. Analyzing problems on multiple scales is 
common in many fields such as CFE^and image processing® . 22 Full multigrid solution methods have even 
been applied to seismic data deconvolution. 23 ^ Currently, only the simple case is explored, whereby coarse 
scales are used to initialize fine scales. 

The construction of scales for deconvolution analysis requires some care. If the coarse grids have very 
few points per beamwidth, the deconvolution results may be unreliable. This can be due to both the limited 
resolution of the array response’s main lobe, as well as the discrete sampling pattern missing side lobes. In 
an attempt to mitigate this coarse-grid limitation, formulations for an average wavespace pressure field and 


average array response are derived and used instead of those given by Eqs. (12) and (14). 
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For the moment, consider a one-dimensional wavespace transform of the pressure field measured by a 
linear array, 


1 M 

p( kx ) = E Pme 


—jk x Xr, 


(26) 


m= 1 


For a grid spacing A fc, the average value of the wavespace transform at grid point k XiU is 


_ E,u+^fe/2 i M 


L 


u +A k /2 e -jk x 


m= 1 
M 


k XjU —Ak/2 


Ak 


dk x . 


(27) 


— ]\f ^ {kx, U i Afe, X m ) 


When = 0, 4> = 1. When 7 ^ 0, 

4^ {J^X^Ul Afc, x m ) = 


m=l 


1 


o 3 k x x n 


-jx m Ak 


fc*,u+Afc/2 


k x : u Ak 1 2 


(28) 


-jx m Ak l 


0 -jXm(k XjU -\-Ak/2) —jx rn (k XjU —Ak/2) 


Using the definition from Eq. (28), the average 4-dimensional wavespace cross-spectrum is constructed as 

_ 1 M M 

P (^,U5 kr^ u ^ k x , u , ky, u ) = 2 ^ ^ ^ ] P [PnPm] X (29) 

n=l m=l 

^ Afc, «^ n ) 4/ (^ )U , Afc, 7y n ) 4/ {k x , U i Afc, x m ) 4/ (fcy ?n , Afc, 2 / 777 ,) 5 

where the negative sign on the conjugate-dimension coordinates in 4/ accounts for the conjugate phase of 
the integration. Similarly, the average cross-power sampling function is given by 


P (^,W) ^?7,U5 ^X,U1 fcj/, 


1 


M M 
n=l m=l 

4^ (fc^,U 5 Afc, ^ n ) 4/ ( ky^u , Afc, Tj n) 4/ {k x ,w> Afc, X m ) 4/ ( fcy,n, Afc, 2 / 777 ,) • 


(30) 


As with the non- averaged formulations, the M 2 denominator term is replaced with M 2 — M when diagonal 
removal is applied and m = n terms neglected. As a final correction, both P and S' are normalized by 


S (0,0, 0,0), since this value is no longer unity except for very fine grids (where Eqs. (29) and (30) collapse 
to Eqs. ( [I2| ) and (14), respectively). As the correction is effectively applied to both sides of Eq. |1J, equality 
of the equivalent system of equations is maintained. 

With modified equations defined, a means of transitioning from more coarse grids to more fine grids must 
now be determined. If some knowledge of the source field Q is available a-priori, a model-based interpolation 
scheme would be preferable. However, in constructing a general deconvolution method, such knowledge 
cannot be assumed. For example, in the case of discrete plane waves, the source field is discontinuous. To 
mitigate the influence of discontinuities on blind interpolation, a linear interpolation scheme is considered. 
To account for potential incoherent source regions in the interpolation scheme, autospectral terms of Q are 
upscaled separately from cross-spectral terms. This interpolation dependency is shown schematically for 
a two-dimensional conjugate-power case in Figure [l] where a 3x3 grid is upscaled to a 5x5 grid. In this 
representation, the wavespace autospectrum exists on the grid diagonal, with cross-spectral terms elsewhere. 
Note that conjugate symmetry exists between the upper and lower triangular portions of the grid. 


V. Plane Wave Results 

Results for simulated plane wave fields, which are identically shift-invariant, are presented. Fields of 
varying complexity are considered. All simulations use the microphone layout of the outer UFAFF array 1 — 
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(a) Pre-interpolation grid (b) Interpolation dependencies 

Figure 1. Linear grid interpolation for 2-D problem upscaling. 


mentioned in Section III. Details of the array can be found in the reference, but in brief the array is a 0.74 
m aperture, 5-arm log spiral based on the methodology of Underbrink^ consisting of 40 total microphones. 
The array’s average autospectrum sampling function S (fc x , k y ,k x ,k y ) is shown for two different resolutions 
in Figure |2j The 3-dB beamwidth of the autospectrum of the sampling function is 1.47 m _1 . The coarse 
grid with n s = 17 points thus has 3 points per beamwidth, while the fine grid with n s = 33 points has 5 
points per beamwidth. Diagonal removal is applied for all cases. For array response plots, the color scale is 
in dB referenced to the array response center value of unity. For plots of pressure and source spectra, the 
color scale is in dB ref 20 f. iPa . As pressure autospectra can contain negative values with diagonal removal, 
these values are set to — oo dB for plotting. 



(a) n s = 17 points (b) n s = 33 points 

Figure 2. Autospectrum of average sampling function S for two different resolutions. 


Isolated Plane Wave 

The first case considered is a single, broadside plane wave, corresponding to a point source in wavespace with 
k x = 0 m -1 and k y = 0 m -1 . The plane wave is simulated with an amplitude of 100 dB. The simulation 
parameters use the more coarse sampling grid from Figure |2a| corresponding to = 9 points and n s = 17 
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points, with k max = 2.92 m . The average wavespace auto-spectrum P ( k x , k y ,k x ,k y ) for this case is plotted 
in Figure [3a] As expected for a centered point source, it appears as a scaled version of the average sampling 


function. The average cross-spectrum magnitude 


P (0, 0, fc x , ky 


is shown in Figure 


After 10,000 iterations, single-precision convergence as defined by Eq. ( [19] ) occurs, requiring 100 seconds 
of runtime on a 2.2 GHz Core i7 MacBook Pro for code utilizing the FFTW engine present in MATLAB 


3b 


2011b. The resultant auto- and cross-spectra, Q (fc x , k y , k x , k y ) and Q (0, 0, k x , k y ) 


, are shown in Figures 


and [3d] As expected, only one source is observed in the spectrum and no cross-powers are evident. 


3c 




k (m 1 ) 

X x 7 

(c) Source autospectrum Q (k x , k y ,k x ,k y ) 



(d) Source cross-spectral magnitude \Q (0,0, k x , k y ) 


Figure 3. Broadside plane wave pressure and source spectra. 


Coherent Plane Wave Pair 


A two source configuration is now analyzed. Here, a field of two plane waves arriving from opposite directions 
is simulated. The waves are defined to correspond to k x = ±1.46 m -1 , k y = 0 m -1 . These values are 
selected to co- locate with grid-point centers. The plane waves are defined as being perfectly coherent, with 
the k x = 1.46 m -1 wave having a 7r/2 radian phase lag behind the k x = —1.46 m -1 wave. Wave magnitudes 
are set to 100 dB. The average wavespace auto- and cross-spectra are respectively plotted in Figures |4a| and 
4b| The coordinates of the reference point for the cross-spectrum are selected to show the relationship of 


the first wave to the second. 
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Convergence is reached after 150,000 iterations, showing a dramatic increase over the single point source 
case. Computation time scales similarly. The resultant source auto- and cross-spectra are shown in Figures 
[4c]and[4d| The source plots again behave as expected for perfectly-coherent sources, where the cross-spectral 
magnitude matches the autospectral magnitudes for equal-power sources. The phase angle between the two 
sources ZQ (—1.46, 0, 1.46, 0) matches 7r/2 to the third decimal place with a value of 1.5709 radians, or a 
relative error of 8.7 x 10 -5 . 




(c) Source autospectrum Q (k x , k y ,k x ,k y ) 



(d) Source cross-spectral magnitude \Q (—1.46, 0, k x , k y ) 


Figure 4. Coherent plane wave pair pressure and source spectra, = 9. 


Multiscale 

This solution is upscaled to an = 13 point grid (4 points per beamwidth) for the same wavespace extent 
kmax = 2.92 m -1 and benchmarked. 1.5 million iterations are required for convergence at this scale. The 
upscaling method is attempted again with a reduction parameter of ( Ak nk= is ) / (A& n/c =g) . Since Q is a 
4-D space of 2-D power spectra, this is equivalent to scaling power-spectral densities. This method fares no 
better with regard to numeric convergence, but does show correct results visually in fewer iterations. 

Conversely, when the solver is started from a Q = 0 condition for the same grid, only 500,000 iterations 
are required for convergence. This type of behavior is observed for most of the discrete, ideal plane wave cases 
considered. Based on observation of intermediate solution states, upscaling appears to bring the solution 
close to the correct value far more quickly than a zero-condition start. However, the solver takes a very long 
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time to completely remove interpolation artifacts, so numerical convergence requires more iterations than 
the zero-start condition. Upscaling by the scaled solution reduces the influence of the interpolation artifacts, 
but shifts the initial guess of the correct bin further from the solution. This is sensible, as interpolating 
these single-wavenumber spatial tones to a more refined grid can be considered analogous to interpolating a 
temporal tone to a more refined spectrum with power-spectral density scaling. The level shift from the scaling 
causes a significant offset in the tone level. Upscaling is thus not considered for the remaining ideal plane 
wave simulations, as these are all defined as discrete spatial tones. However, power-spectral density-based 
upscaling is revisited in Section VI, where distributed field behavior is evident. 

The 500,000 iteration zero-start case on the rik = 13 grid required just under 3 hours of run-time. 
The average pressure auto- and cross-spectra are respectively plotted in Figures |5a| and |5b| The resultant 
source auto- and cross- spectra are shown in Figures [5c] and |5d| The phase angle between the two sources 
ZQ (—1.46, 0, 1.46, 0) shows slightly worse agreement with the true input, as it has a value of 1.5697 radians, 
or a relative error of 7.0 x 10 -4 . 




k (m' 1 ) 

X V ' 

(c) Source autospectrum Q (k x , k y ,k x ,k y ) 



(d) Source cross-spectral magnitude \Q (—1.46, 0, k x , k y ) 


Figure 5. Coherent plane wave pair pressure and source spectra, rik = 13. 


Mixed-Coherence Isolated Plane Waves 

On the rik = 9 point grid, this two source configuration is further complicated with the addition of two 
incoherent sources to the field. The new plane waves again arrive from opposite directions, this time defined 
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at k x = 0 m , k y = ±1.46 m . The average auto- and cross-spectra are plotted in Figures 6a and 6b 
The resultant source auto- and cross-spectra are shown in Figures [6c] and [6d] Numeric convergence requires 
500,000 iterations. As shown, the method properly separates the coherent and incoherent sources in the 
field. The phase angle ZQ (—1.46, 0, 1.46, 0) is again matched to the third decimal place with a value of 
1.5706 radians, or a relative error of 1.2 x 10 -4 . 




(c) Source autospectrum Q (k x , k y ,k x ,k y ) 



(d) Source cross-spectral magnitude \Q (—1.46, 0, k x , k y ) 


Figure 6. Pressure and source spectra for pair of coherent and pair of incoherent sources. 


Line Source with Decaying Coherence 


The final ideal plane wave case considered consists of a line of sources located along the k x axis, at values of 
k x = —1.46, —0.73, 0, 0.73, 1.46 m -1 and k y = 0 m -1 . As before, all of the source auto-powers are set to 100 
dB. However, the sources are defined to have decaying coherence relative to the center of the line. Defining 
a coherence-squared function as 


T (&£ 5 ky , k x , ky') 


Q {k^t ky , k x , ky ) 


Q (^£, ^775 ^77) Q (Mr 5 kyi •> ky) 


(31) 


the values referenced to the center are set to y 2 (0, 0, ±0.73, 0) = 0.5 and y 2 (0, 0, ±1.46, 0) = 0.25. The phase 
angles of the coherent powers relative to the center source are set to ZQ (0, 0, ±0.73, 0) = ±7r/4 radians and 


14 ofl23l 


American Institute of Aeronautics and Astronautics 


ZQ (0, 0, ±1.46, 0) = =F7r/2 radians. 

The average pressure autospectrum is plotted in Figure [7a] The average pressure cross-spectral magnitude 
is plotted in Figure [7b] A strong interaction between the sources is present, as the peak value of the pressure 
field is 6 dB higher than any of the input sources. This interaction is evident in both the auto- and cross- 
spectra. 

The algorithm is run for 5 million iterations and numerical convergence is not achieved, although residual 
levels are close to the convergence value. There is clearly a very strong convergence dependence on source 
field complexity. However, the results after 5 million iterations, plotted for the source field autospectrum 
in Figure |7c| and cross-spectral magnitude in Figure |7d] show the expected source behavior, aside from one 
spurious source 20 dB below peak in the cross-spectrum. All of the autospectral values are at 100 dB, and 
the cross- spectral values decay from 100 dB at k x = 0 m -1 to 97 dB at k x = ±1.46 m -1 . 





k (m 1 ) k (m -1 ) 

x v 7 x x 7 

(c) Source autospectrum Q ( k x , k y ,k x ,k y ) (d) Source cross-spectral magnitude | Q (0, 0, k x ,k y ) 

Figure 7. Pressure and source spectra for a line of partially-coherent sources. 


Figure [8] shows the behavior of the line of sources in more detail. Figure [8a] plots the cross-spectral 


magnitude 


Q(%, 0,^,0) 


as all of the non-zero elements of the true solution exist on this plane. As 
expected, the magnitude decays away from the autospectrum, which exists on the diagonal of the grid. 
Figure [8b] takes the k^ = 0 cut of this plane and compares it to the true solution. As shown on a decibel 
scale, the results are close. Figure [8c] then shows the cross-spectrum normalized to the coherence-squared and 


compares to the input values. Finally, Figure 8d shows the resultant phase angles ZQ (0, 0, k x , 0) compared 
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to the input angles. From this simulation, it is apparent that true numerical convergence is not necessary 
for reasonable results. However, not shown are the partial results from the first 2 million iterations, where 
spectral levels were still over 1 dB off, and phase angle errors were greater than 5%. 
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Figure 8. Evaluation of line source results and comparison to inputs. 


VI. Curved- Wave Beamforming 

Results from Section V show that this deconvolution formulation is functional, if expensive, for ideal 
discrete plane waves of varying coherence. However, under many circumstances an aeroacoustic array will 
observe an acoustic field which consists of many curved wave fronts. While a true curved wave front cannot 
be fully represented by a select few frequencies in wavespace, it may be possible to find a “best-fit” wavespace 
spectrum to a given curved wave field. 

Simulated Results 

An initial simulation consists of a point source, centered over the array, located 1.5 m away from the array 
center, for a source coordinate of (0,0, 1.5). This location is selected as approximately twice the array aperture 
to generate a moderately- curved wavefront as observed by the acoustic array. The source is scaled such that 
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its acoustic field has a level of 90.5 dB at the array center, and 90.2 dB at the outer microphone ring. As 
the acoustic field of this point source will appear as the superposition of many plane waves, the wavespace 
spectrum must be constructed to capture the entire acoustic radiation circle. If k max from the previous 
sections is maintained and set to the radius of the radiation circle, computed as k max = //cop^this sets the 

source frequency at 1 kHz. The average spectrum P is constructed on three grid scales of rik = 9, n k = 13 
and rik = 17. For constant & maa ,, this gives successively smaller A k values, corresponding to increasingly 
large effective array apertures. Numeric convergence is not expected for these cases based on previous 
results, so the algorithm is run until residuals and error no longer decrease. Figure |9a| shows the average 

shows the solution autospectrum after 500,000 


9b 


pressure autospectrum P for the rik = 9 grid, while Figure]! 
iterations. This solution is then upscaled to the n k = 13 grid. The n & = 13 average pressure autospectrum 
and solution autospectrum after 50,000 iterations are shown in Figures 9c and 9d This solution is then 


upscaled to the *= 17 grid, with average pressure spectrum and solution after 20,000 iterations shown in 
Figures [9e] and [9f[ 

The plotted wavespace results appear as a cone of arrival directions centered about the broadside direction, 
as expected. For this simulation, the array occupies a 28° solid angle of the source’s radiation pattern. While 
explicit definition of the true solution for a spatially-varying, array-sampled, curved wavefront in wavespace 
may be difficult, it would not be unreasonable to expect that the source field have the majority of its power 
occupy a direction-of-arrival cone corresponding to this solid angle. For this case, the energy should reside 
approximately within a circle of radius k = 0.7 m -1 , which is seen for all three rik values. Quantitatively, 
the total summed power within the source autospectra is 89.5 dB for n & = 9, 88.9 dB for n k = 13, and 88.4 
dB for rik = 17. These values do not match the expected levels for the acoustic power observed by the array. 
However, if the total source field is summed over all four wavespace dimensions, the total coherent power 
sum can be computed. This total summed power within the wavespace spectrum is 90.5 dB for all three 
grid densities, matching the true acoustic field power at the array center. The rik = 17 grid is restarted with 
a zero-source field to evaluate the effects of upscaling. This case is run for 20,000 iterations. For brevity, 
no plots are shown. The total summed power of the source field for this case is still 90.5 dB. However, the 
summed autospectral power is only 80.9 dB. While there is no true solution to compare and say which result 
is closer, the error as defined in Equation © is two orders of magnitude smaller for the upscaled case than 
for the zero-start case. 


Experimental Results: Trailing Edge Noise 

The final case considered is taken from an existing trailing edge noise data set collected in UFAFF. 26 Data 
are acquired for a 0.3048 m chord, 0.74 m span NACA 0012 airfoil, in this case at a Mach number of 0.17. 
Array data are acquired for an array with the same layout as that used in all previous simulations. A laser 
pulse calibration technique is applied to the data CSM to calibrate the individual microphone channels and, 
to first order, correct for shear layer refraction effects? 6 Additional details of the experiment can be found 
in the reference. A photograph of the installation is shown in Figure [To] while a legend of the illustrations 
used in the baseline beam map at z = 1.13 m is shown in Figure El Here, flow is from right to left. The 
airfoil trailing edge center span is located at (0,0,1.13) m from the array center. A conventional beam map of 
the array data is shown in Figure [12] Positive x in the beam map points in the upstream direction. Positive 
y points upward vertically, while positive z points out of the page. The airfoil trailing edge appears as the 
dominant noise source. 

The wavespace deconvolution procedure is applied to the 2 kHz narrowband CSM computed from the 
experiment. This is the only frequency bin considered for the following analysis. This increase in frequency 
from 1 kHz to 2 kHz requires that, to maintain beamwidth resolution, rik be increased to 33 points to account 
for the increase in acoustic wavenumber and & maa ,. 

The algorithm is run on a succession of grids of rik = 9, rik = 13, rik — 17, rik = 25 and rik = 33. A 
surprising effect is observed for this case. While in previous cases the algorithm required a large number 
of iterations at every stage, here a clean, logarithmic roll-off in the error and residual from Eqs. (19) and 
(20) is observed, and relatively few iterations are required for near-total single-precision convergence. At 
rik = 9, 10,000 iterations are required. However, at both rik = 13 and rik = 17, only 3,000 are necessary. 
At rik = 25 15,000 are required for convergence, and finally 20,000 are needed for rik = 33. This requires 
just over 7 hours of computation time. The validity of this upscaling is evaluated by running the rik = 17 
grid from a zero- start. The autospectrum quickly attains a very similar shape to that seen in the upscaling 
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(c) P {kxt ky, k X: ky), rife — 13 


(d) Q (k x , k y , k x , k y ), n k = 13 



Figure 9. Pressure and source autospectra for a near-array point source. 
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Figure 10. NACA 0012 trailing edge noise measurement installation. 



Figure 11. Facility boundary lines as drawn in beam map. 
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Figure 12. Baseline beamforming for NACA 0012 experiment, z = 1.13 m. 


process, but levels rise slowly and, after an initial drop, the error and residual behavior is worse. There are 
many possible reasons for this. The effect may indicate that the interpolation scheme is reasonable for this 
case, or that this particular wave field lends itself readily to representation by a summation of plane waves. 
Alternately, the sampling grid may have better conditioning, or there could be some other numerical effect 
such as floating point precision which did not matter for smaller grids, but comes into play here. Further 
investigation is necessary. 

The resultant average and source autospectra for the = 9, = 17 and = 33 grids are shown in 

Figure [13] Broadside waves at k x = 0 are clearly present, as expected for an airfoil trailing edge acoustic 
source located at x = 0 m. Additionally, waves arriving from downstream, located near k x = —4 m -1 , are 
likely due to the facility diffuser. Additional low-level sources appear outside of the radiation circle. Based 
on the average pressure spectra, these may not be sensible. Additional comparison of upscaled and zero-start 
solutions is likely required. The total sum power of the wavespace spectrum for the trailing edge region of 
the source field for = 33 is 45.7 dB. The sum power of the autospectrum in this region and at this scale is 
37.6 dB. The two- microphone method results from the reference, concluded to be the most accurate estimate 
of the trailing edge noise for this frequency, showed a trailing edge noise level of 48 dB. 

Proposed Extension to Nearfield Beamforming 

Initial results for curved waves indicate that approximate source region powers may be computed by sum- 
ming the appropriate subset of wavespace auto- and cross-powers for given directions-of- arrival of interest. 
However, interpretation of these results is, at the moment, difficult and potentially non-intuitive depending 
on source configuration. However, with a full cross-power spectrum in wavespace, a synthetic CSM may be 
constructed by inverse-transforming the spectrum Q. This CSM can be used with conventional beamform- 
ing in an attempt to further isolate sources, and possibly with cross-beamformin^ to determine regional 
coherence. Initial tests with this process show qualitative success, but are far from quantitative reliability. 
Further investigation of this is left to future work. 

VII. Summary and Future Work 

A coherent deconvolution technique is presented for a class of shift- invariant problems constructed in 
wavespace. The deconvolution method is based on an iterative solver which enforces appropriate constraints 
after every iteration. The solver iteratively evaluates a forward convolution utilizing a four-dimensional FFT, 
providing improved scaling with problem size when compared to matrix methods. The technique is applied 
to ideal plane wave cases successfully, and properly recovers both magnitude and phase of cross-spectral 
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Figure 13. Pressure and source autospectra for the NACA 0012 experiment. 
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relationships between distinct waves of arbitrary coherence. Acceleration through an improved calculation 
of a stable relaxation parameter greatly improves the behavior of the algorithm. Upscaling for an initial 
guess on a given grid performs poorly for discrete plane waves when combined with linear interpolation. 
However, it appears to provide an improved initial guess for true distributed noise sources. 

A large amount of work remains to properly evaluate this method. First and foremost, more curved 
wave simulations must be conducted with a variety of multi-source configurations. The reconstruction of 
synthetic CSM s must be pursued. The upscaling procedure must be further investigated with regard to the 
convergence behavior of the trailing edge noise data set. Additionally, the effect of the point source shear 
layer correction on the data is not considered here, and may affect solution behavior. Experimental data 
which do not require such correction should be analyzed in any future study. The effect of floating-point 
precision on larger grids with real data should be considered. 

Acceleration must continue to be pursued. While FFTs show favorable scaling with advances in compu- 
tational technology, alternative means of improving solution times must still be sought. Multigrid methods 
may prove to be the most promising for distributed noise fields, but four-dimensional interpolation of cross- 
spectra must be evaluated in more detail, given the observed behavior with discrete plane waves. Careless 
application of multigrid or even simple upscaling may prove problematic if both distributed noise sources and 
discrete plane waves are present in a given acoustic field. Alternatively, other algorithms may be explored. 
The iterative technique considered here is one of many ways to pursue what amounts to a constrained opti- 
mization problem. Finding a technique which does not require tens or hundreds of thousands of iterations 
to reach a low-residual solution would be of huge benefit to this deconvolution method. Finally, while this 
wavespace deconvolution algorithm circumvents the full computation and storage of the coefficient matrix A, 
the problem still shows poor scaling for high frequencies, especially with regards to memory consumption. 
An alternative method of performing the forward convolution which avoids the storage and computation 
associated with four-dimensional zero-padding could show large improvements in both computation time 
and memory requirements. 
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